# Result-Plots from Signorino DGP

rm(list = ls())
setwd("/Users/lleemann/Desktop/Replication Files/Simulation/Alternative DGP")

load("ResultSimulation_-09_250.RData")
objects()

param.struc__09 <-  stuff$param.struc
param.selec__09 <-  stuff$param.selec
param.selStruc__09 <-  stuff$param.selStruc

bias.struc__09 <-  stuff$bias.struc
bias.selec__09 <-  stuff$bias.selec
bias.selStruc__09 <-  stuff$bias.selStruc

load("ResultSimulation_-05_250.RData")
objects()

param.struc__05 <-  stuff$param.struc
param.selec__05 <-  stuff$param.selec
param.selStruc__05 <-  stuff$param.selStruc

bias.struc__05 <-  stuff$bias.struc
bias.selec__05 <-  stuff$bias.selec
bias.selStruc__05 <-  stuff$bias.selStruc

load("ResultSimulation_00_250.RData")
objects()

param.struc_00 <-  stuff$param.struc
param.selec_00 <-  stuff$param.selec
param.selStruc_00 <-  stuff$param.selStruc

bias.struc_00 <-  stuff$bias.struc
bias.selec_00 <-  stuff$bias.selec
bias.selStruc_00 <-  stuff$bias.selStruc



load("ResultSimulation_+05_250.RData")
objects()

param.struc_05 <-  stuff$param.struc
param.selec_05 <-  stuff$param.selec
param.selStruc_05 <-  stuff$param.selStruc

bias.struc_05 <-  stuff$bias.struc
bias.selec_05 <-  stuff$bias.selec
bias.selStruc_05 <-  stuff$bias.selStruc



load("ResultSimulation_+09_250.RData")
objects()

param.struc_09 <-  stuff$param.struc
param.selec_09 <-  stuff$param.selec
param.selStruc_09 <-  stuff$param.selStruc

bias.struc_09 <-  stuff$bias.struc
bias.selec_09 <-  stuff$bias.selec
bias.selStruc_09 <-  stuff$bias.selStruc





## 1st plot
k <- 6
bias.struc <- c(mean(param.struc__09[,k]-1), mean(param.struc__05[,k]-1), mean(param.struc_00[,k]-1), mean(param.struc_05[,k]-1), mean(param.struc_09[,k]-1))
bias.selec <- c(mean(param.selec__09[,k-1]-1), mean(param.selec__05[,k-1]-1), mean(param.selec_00[,k-1]-1), mean(param.selec_05[,k-1]-1), mean(param.selec_09[,k-1]-1))
bias.selStruc <- c(mean(param.selStruc__09[,k]-1), mean(param.selStruc__05[,k]-1), mean(param.selStruc_00[,k]-1), mean(param.selStruc_05[,k]-1), mean(param.selStruc_09[,k]-1))

corr <- c(-.9,-.5,0,.5,.9)

#par(family="CMU Serif", mfrow=c(1,3), mar=c(5.1, 2.1, 2.1, 0.5), oma=c(0, 2, 1, 0))
par(family="CMU Serif", mfrow=c(2,3), mar=c(4.1, 4.1, 2.1, 0.5), oma=c(0, 2, 1, 0))
plot(corr,bias.struc,xlim=c(-1,1), ylim=c(-0.12,0.45), pch=16, col="grey60", cex=2.6, ylab="", xlab="", main=expression(paste("Estimate of ", beta[21])))
abline(h=0, lty=2)
points(corr,bias.selec,pch=15, col="grey80", cex=2.6)
points(corr,bias.selStruc,pch=17, col="black", cex=2.6)

# yellow
for (j in 2:6){
		segments(x0=corr[j-1], y0=bias.struc[j-1], x=corr[j], y=bias.struc[j], lty=2, col="grey60", lwd=2)
	}
# green
for (j in 2:6){
		segments(x0=corr[j-1], y0=bias.selec[j-1], x=corr[j], y=bias.selec[j], lty=2, col="grey80", lwd=2)
	}
# purple
for (j in 2:6){
		segments(x0=corr[j-1], y0=bias.selStruc[j-1], x=corr[j], y=bias.selStruc[j], lty=2, col="black", lwd=2)
	}

legend(.4,.31, legend=c("Strat", "Heck", "SSa"), col=c(col="grey60","grey80","black"), cex=1.2, pch=c(16,15,17), bty="n")


## 2nd plot
k <- 7
bias.struc <- c(mean(param.struc__09[,k]-1), mean(param.struc__05[,k]-1), mean(param.struc_00[,k]-1), mean(param.struc_05[,k]-1), mean(param.struc_09[,k]-1))
bias.selec <- c(mean(param.selec__09[,k-1]-1), mean(param.selec__05[,k-1]-1), mean(param.selec_00[,k-1]-1), mean(param.selec_05[,k-1]-1), mean(param.selec_09[,k-1]-1))
bias.selStruc <- c(mean(param.selStruc__09[,k]-1), mean(param.selStruc__05[,k]-1), mean(param.selStruc_00[,k]-1), mean(param.selStruc_05[,k]-1), mean(param.selStruc_09[,k]-1))


plot(corr,bias.struc,xlim=c(-1,1), ylim=c(-0.12,0.45), pch=16, col="grey60", cex=2.6, ylab="", xlab="", main=expression(paste("Estimate of ", beta[22])))
abline(h=0, lty=2)
points(corr,bias.selec,pch=15, col="grey80", cex=2.6)
points(corr,bias.selStruc,pch=17, col="black", cex=2.6)

# yellow
for (j in 2:6){
		segments(x0=corr[j-1], y0=bias.struc[j-1], x=corr[j], y=bias.struc[j], lty=2, col="grey60", lwd=2)
	}
# green
for (j in 2:6){
		segments(x0=corr[j-1], y0=bias.selec[j-1], x=corr[j], y=bias.selec[j], lty=2, col="grey80", lwd=2)
	}
# purple
for (j in 2:6){
		segments(x0=corr[j-1], y0=bias.selStruc[j-1], x=corr[j], y=bias.selStruc[j], lty=2, col="black", lwd=2)
	}

legend(.4,.31, legend=c("Strat", "Heck", "SSa"), col=c(col="grey60","grey80","black"), cex=1.2, pch=c(16,15,17), bty="n")


## 3rd plot --b12
k <- 2
bias.struc <- c(mean(param.struc__09[,k]-1), mean(param.struc__05[,k]-1), mean(param.struc_00[,k]-1), mean(param.struc_05[,k]-1), mean(param.struc_09[,k]-1))
bias.selec <- c(mean(param.selec__09[,k+1]-1), mean(param.selec__05[,k+1]-1), mean(param.selec_00[,k+1]-1), mean(param.selec_05[,k+1]-1), mean(param.selec_09[,k+1]-1))
bias.selStruc <- c(mean(param.selStruc__09[,k]-1), mean(param.selStruc__05[,k]-1), mean(param.selStruc_00[,k]-1), mean(param.selStruc_05[,k]-1), mean(param.selStruc_09[,k]-1))


plot(corr,bias.struc,xlim=c(-1,1), ylim=c(-0.62,0.25), pch=16, col="grey60", cex=2.6, ylab="", xlab="", main=expression(paste("Estimate of ", beta[12])))
abline(h=0, lty=2)
points(corr,bias.selec,pch=15, col="grey80", cex=2.6)
points(corr,bias.selStruc,pch=17, col="black", cex=2.6)

# yellow
for (j in 2:6){
		segments(x0=corr[j-1], y0=bias.struc[j-1], x=corr[j], y=bias.struc[j], lty=2, col="grey60", lwd=2)
	}
# green
for (j in 2:6){
		segments(x0=corr[j-1], y0=bias.selec[j-1], x=corr[j], y=bias.selec[j], lty=2, col="grey80", lwd=2)
	}
# purple
for (j in 2:6){
		segments(x0=corr[j-1], y0=bias.selStruc[j-1], x=corr[j], y=bias.selStruc[j], lty=2, col="black", lwd=2)
	}

legend(.4,-.065, legend=c("Strat", "Heck", "SSa"), col=c(col="grey60","grey80","black"), cex=1.2, pch=c(16,15,17), bty="n")





mtext(c("Bias","Root Mean Squared Error"),outer=TRUE,line=.1, side=2, at=c(.78,2))


####	####	####	####	####	####	####	####	####	####	####	####

####	####	####	####	####	####	####	####	####	####	####	####

####	####	####	####	####	####	####	####	####	####	####	####

####	####	####	####	####	####	####	####	####	####	####	####	

### 			2nd ROW  of plots (RMSE)

####	####	####	####	####	####	####	####	####	####	####	####	

####	####	####	####	####	####	####	####	####	####	####	####

####	####	####	####	####	####	####	####	####	####	####	####

####	####	####	####	####	####	####	####	####	####	####	####

## 1st plot
k <- 6
rmse.struc <- c(sqrt(mean((param.struc__09[,k]-1)^2)), sqrt(mean((param.struc__05[,k]-1)^2)), sqrt(mean((param.struc_00[,k]-1)^2)), sqrt(mean((param.struc_05[,k]-1)^2)), sqrt(mean((param.struc_09[,k]-1)^2)))
rmse.selec <- c(sqrt(mean((param.selec__09[,k-1]-1)^2)), sqrt(mean((param.selec__05[,k-1]-1)^2)), sqrt(mean((param.selec_00[,k-1]-1)^2)), sqrt(mean((param.selec_05[,k-1]-1)^2)), sqrt(mean((param.selec_09[,k-1]-1)^2)))
rmse.selStruc <- c(sqrt(mean((param.selStruc__09[,k]-1)^2)), sqrt(mean((param.selStruc__05[,k]-1)^2)), sqrt(mean((param.selStruc_00[,k]-1)^2)), sqrt(mean((param.selStruc_05[,k]-1)^2)), sqrt(mean((param.selStruc_09[,k]-1)^2)))

corr <- c(-.9,-.5,0,.5,.9)

plot(corr,rmse.struc,xlim=c(-1,1), ylim=c(-0.02,0.45), pch=16, col="grey60", cex=2.6, ylab="", xlab="Correlation of Errors")
abline(h=0, lty=2)
points(corr,rmse.selec,pch=15, col="grey80", cex=2.6)
points(corr,rmse.selStruc,pch=17, col="black", cex=2.6)

# yellow
for (j in 2:6){
		segments(x0=corr[j-1], y0=rmse.struc[j-1], x=corr[j], y=rmse.struc[j], lty=2, col="grey60", lwd=2)
	}
# green
for (j in 2:6){
		segments(x0=corr[j-1], y0=rmse.selec[j-1], x=corr[j], y=rmse.selec[j], lty=2, col="grey80", lwd=2)
	}
# purple
for (j in 2:6){
		segments(x0=corr[j-1], y0=rmse.selStruc[j-1], x=corr[j], y=rmse.selStruc[j], lty=2, col="black", lwd=2)
	}

legend(.4,.365, legend=c("Strat", "Heck", "SSa"), col=c(col="grey60","grey80","black"), cex=1.2, pch=c(16,15,17), bty="n")

## 2nd plot, 2nd row
k <- 7
rmse.struc <- c(sqrt(mean((param.struc__09[,k]-1)^2)), sqrt(mean((param.struc__05[,k]-1)^2)), sqrt(mean((param.struc_00[,k]-1)^2)), sqrt(mean((param.struc_05[,k]-1)^2)), sqrt(mean((param.struc_09[,k]-1)^2)))
rmse.selec <- c(sqrt(mean((param.selec__09[,k-1]-1)^2)), sqrt(mean((param.selec__05[,k-1]-1)^2)), sqrt(mean((param.selec_00[,k-1]-1)^2)), sqrt(mean((param.selec_05[,k-1]-1)^2)), sqrt(mean((param.selec_09[,k-1]-1)^2)))
rmse.selStruc <- c(sqrt(mean((param.selStruc__09[,k]-1)^2)), sqrt(mean((param.selStruc__05[,k]-1)^2)), sqrt(mean((param.selStruc_00[,k]-1)^2)), sqrt(mean((param.selStruc_05[,k]-1)^2)), sqrt(mean((param.selStruc_09[,k]-1)^2)))

corr <- c(-.9,-.5,0,.5,.9)

#par(family="CMU Serif", mfrow=c(2,3), mar=c(5.1, 4.1, 1.1, 0.5))
#par(family="CMU Serif")
plot(corr,rmse.struc,xlim=c(-1,1), ylim=c(-0.02,0.45), pch=16, col="grey60", cex=2.6, ylab="", xlab="Correlation of Errors")
abline(h=0, lty=2)
points(corr,rmse.selec,pch=15, col="grey80", cex=2.6)
points(corr,rmse.selStruc,pch=17, col="black", cex=2.6)

# yellow
for (j in 2:6){
		segments(x0=corr[j-1], y0=rmse.struc[j-1], x=corr[j], y=rmse.struc[j], lty=2, col="grey60", lwd=2)
	}
# green
for (j in 2:6){
		segments(x0=corr[j-1], y0=rmse.selec[j-1], x=corr[j], y=rmse.selec[j], lty=2, col="grey80", lwd=2)
	}
# purple
for (j in 2:6){
		segments(x0=corr[j-1], y0=rmse.selStruc[j-1], x=corr[j], y=rmse.selStruc[j], lty=2, col="black", lwd=2)
	}

legend(.4,.365, legend=c("Strat", "Heck", "SSa"), col=c(col="grey60","grey80","black"), cex=1.2, pch=c(16,15,17), bty="n")


## 3rd plot, 2nd row
k <- 2
rmse.struc <- c(sqrt(mean((param.struc__09[,k]-1)^2)), sqrt(mean((param.struc__05[,k]-1)^2)), sqrt(mean((param.struc_00[,k]-1)^2)), sqrt(mean((param.struc_05[,k]-1)^2)), sqrt(mean((param.struc_09[,k]-1)^2)))
rmse.selec <- c(sqrt(mean((param.selec__09[,k+1]-1)^2)), sqrt(mean((param.selec__05[,k+1]-1)^2)), sqrt(mean((param.selec_00[,k+1]-1)^2)), sqrt(mean((param.selec_05[,k+1]-1)^2)), sqrt(mean((param.selec_09[,k+1]-1)^2)))
rmse.selStruc <- c(sqrt(mean((param.selStruc__09[,k]-1)^2)), sqrt(mean((param.selStruc__05[,k]-1)^2)), sqrt(mean((param.selStruc_00[,k]-1)^2)), sqrt(mean((param.selStruc_05[,k]-1)^2)), sqrt(mean((param.selStruc_09[,k]-1)^2)))

corr <- c(-.9,-.5,0,.5,.9)

#par(family="CMU Serif", mfrow=c(2,3), mar=c(5.1, 4.1, 1.1, 0.5))
#par(family="CMU Serif")
plot(corr,rmse.struc,xlim=c(-1,1), ylim=c(-0.02,0.6), pch=16, col="grey60", cex=2.6, ylab="", xlab="Correlation of Errors")
abline(h=0, lty=2)
points(corr,rmse.selec,pch=15, col="grey80", cex=2.6)
points(corr,rmse.selStruc,pch=17, col="black", cex=2.6)

# yellow
for (j in 2:6){
		segments(x0=corr[j-1], y0=rmse.struc[j-1], x=corr[j], y=rmse.struc[j], lty=2, col="grey60", lwd=2)
	}
# green
for (j in 2:6){
		segments(x0=corr[j-1], y0=rmse.selec[j-1], x=corr[j], y=rmse.selec[j], lty=2, col="grey80", lwd=2)
	}
# purple
for (j in 2:6){
		segments(x0=corr[j-1], y0=rmse.selStruc[j-1], x=corr[j], y=rmse.selStruc[j], lty=2, col="black", lwd=2)
	}


legend(.4,.445, legend=c("Strat", "Heck", "SSa"), col=c(col="grey60","grey80","black"), cex=1.2, pch=c(16,15,17), bty="n")

mtext("Root Mean Squared Error",outer=TRUE,line=.1, side=2, at=c(0.28))





####			####			####			###
